Fit Tweedie models to biomass density of cod, flounder, plaice and dab (juveniles and adults) between 1993-2020 in the Baltic Sea using sdmMTB fit different oxygen and temperature-related covariates. Predict on grid and use those predictions to calculate trends and velocities.
# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(readxl)
library(tidylog)
library(RCurl)
library(devtools)
library(sdmTMB)
library(rMR)
library(kableExtra)
# Source code for plots
source_url("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/R/functions/map-plot.R")
#> Warning: package 'sf' was built under R version 4.0.5
theme_set(theme_plot())
# Read data
d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/catch_clean.csv") %>%
rename(X = x, Y = y) %>%
drop_na(depth) %>%
mutate(depth_ct = depth - mean(depth),
depth_sc = depth_ct / sd(depth),
depth_sq = depth_sc^2,
year_f = as.factor(year),
quarter_f = as.factor(quarter)) %>%
pivot_longer(c(cod_adult, cod_juvenile, dab_adult, dab_juvenile, flounder_adult,
flounder_juvenile, plaice_adult, plaice_juvenile),
names_to = "group", values_to = "density") %>%
separate(group, into = c("species", "life_stage"), remove = FALSE)
#> Rows: 9376 Columns: 23
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (3): country, haul_id, ices_rect
#> dbl (20): year, lat, lon, quarter, month, sub_div, oxy, temp, sal, depth, x,...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> rename: renamed 2 variables (X, Y)
#>
#> drop_na: removed 112 rows (1%), 9,264 rows remaining
#>
#> mutate: new variable 'depth_ct' (double) with 3,707 unique values and 0% NA
#>
#> new variable 'depth_sc' (double) with 3,707 unique values and 0% NA
#>
#> new variable 'depth_sq' (double) with 3,707 unique values and 0% NA
#>
#> new variable 'year_f' (factor) with 28 unique values and 0% NA
#>
#> new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#>
#> pivot_longer: reorganized (cod_adult, cod_juvenile, dab_adult, dab_juvenile, flounder_adult, …) into (group, density) [was 9264x28, now 74112x22]
# Read metabolic parameter estimates and left_join
# mi_pars <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/mi_params.csv") %>%
# dplyr::select(-...1, -temp)
# TODO: for now we'll use plaice parameters for flounder, see "00_estimate_mi_params.Rmd"
# mi_pars2 <- mi_pars
# mi_pars2[4, c(4:5)] <- mi_pars2[1, c(4:5)]
#
# d <- left_join(d, mi_pars2, by = "species")
# Read size csv to calculate the metabolic index
# sizes <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/sizes.csv") %>%
# mutate(group = paste(species, name, sep = "_")) %>%
# dplyr::select(group, B)
#
# d <- left_join(d, sizes, by = "group")
# Read prediction grid
pred_grid <- bind_rows(readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/master/data/clean/pred_grid_(1_2).csv"),
readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/master/data/clean/pred_grid_(2_2).csv"))
#> Rows: 454412 Columns: 12
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): ices_rect
#> dbl (11): X, Y, year, lon, lat, depth, quarter, oxy, temp, sal, sub_div
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> Rows: 454412 Columns: 12
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): ices_rect
#> dbl (11): X, Y, year, lon, lat, depth, quarter, oxy, temp, sal, sub_div
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# # Before we can calculate the pressure based metabolic index, we'll need to convert oxygen from ml/L mg/L and partial pressure
# # Some conversions
# # oxygen: https://www.ices.dk/data/tools/Pages/Unit-conversions.aspx
# # pressure: https://bluerobotics.com/learn/pressure-depth-calculator/
# # First calculate pressure in pascal: p_tot = (r*g*h) + p_atm, then to kgPa and finally atm
# rho = 1030 # density of water
# p_atm = 1.01325 # atmospheric pressure
# g = 9.807 # gravity of earth
#
# # Convert to partial pressure oxygen from ml/L
# d$po2 <- DO.unit.convert(d$oxy*(1/0.7), # d$oxy is in unit ml/L so we first need to convert to mg/L
# DO.units.in = "mg/L",
# DO.units.out = "PP",
# bar.press = ((rho*g*d$depth + p_atm)*10^-5)*0.986923,
# bar.units.in = "atm",
# bar.units.out = "kpa",
# temp.C = d$temp,
# salinity = 10, # TODO: d$sal, # let's fix salinity for easier comparison. It has a minor effect in this range
# salinity.units = "pp.thou")
#
# # Oxygen is ml/L, We want micro mol/L. 1 ml/l = 10^3/22.391 = 44.661 micro mol/l
# d$oxy_si <- (d$oxy * (10^3)) / 22.391
#
# # Calculate metabolic indices for a given mass and the temperature and oxygen in data
# # Line 123 in https://github.com/fate-spatialindicators/SDM_O2/blob/master/code/mi_functions.R
# kb <- 0.000086173324 # Boltzmann's constant
# t_ref <- 15 # arbitrary reference temperature
#
# # Calculate the metabolic index for pressure and concentration
# d <- d %>%
# mutate(inv_temp = (1/(temp + 273.15) - 1/(t_ref + 273.15)),
# phi = A0_po2*(B^n_po2)*po2 * exp((E_po2/kb)*inv_temp),
# phi_c = A0_o2*(B^n_o2)*oxy_si * exp((E_o2/kb)*inv_temp))
# Fit models and save grid-level predictions
data_list <- list()
for(i in unique(d$group)) {
dd <- d %>%
filter(group == i) %>%
droplevels() %>%
drop_na(density)
mesh <- make_mesh(dd, c("X", "Y"), n_knots = 200)
plot(mesh)
m <- sdmTMB(
density ~ 0 + year_f + quarter_f + depth_sc + depth_sq,
data = dd, mesh = mesh, family = tweedie(link = "log"),
spatial = "on", time = "year", spatiotemporal = "IID",
control = sdmTMBcontrol(newton_loops = 1))
sanity(m)
# Plot residuals
res <- residuals(m) %>% as.data.frame()
ggplot(res, aes(sample = V1)) + stat_qq() + stat_qq_line() + theme(aspect.ratio = 1)
ggsave(paste0("figures/supp/qq_density_sdm/sdm_02_", i, ".pdf"),
width = 17, height = 17, units = "cm")
# Save model object
saveRDS(m, paste("output/models/sdm_02_", i, ".rds", sep = ""))
# Predict on grid
pred_grid <- pred_grid %>%
mutate(year_f = as.factor(year),
depth_ct = depth - mean(dd$depth),
depth_sc = depth_ct / sd(dd$depth),
depth_sq = depth_sc^2,
quarter_f = as.factor(quarter))
pred <- predict(m, newdata = pred_grid) %>% mutate(model = i)
data_list[[i]] <- pred
}
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 13 rows (<1%), 9,251 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'year_f' (factor) with 28 unique values and 0% NA
#> new variable 'depth_ct' (double) with 11,144 unique values and 0% NA
#> new variable 'depth_sc' (double) with 11,144 unique values and 0% NA
#> new variable 'depth_sq' (double) with 11,144 unique values and 0% NA
#> new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 13 rows (<1%), 9,251 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 629 rows (7%), 8,635 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 629 rows (7%), 8,635 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 235 rows (3%), 9,029 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 235 rows (3%), 9,029 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 622 rows (7%), 8,642 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 622 rows (7%), 8,642 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
# Save predictions and sims as data frames
df <- dplyr::bind_rows(data_list)
write_csv(df, "output/grid_pred_biomass_02.csv")